# helper function
rowDiff <- function(mat){
  mat[,1] <- mat[,1]*-1
  return(rowSums(mat))
}

# function to get credible interval from samples
ci_bb_sim <- function(s_id, qtile = c(0.05,0.5,0.95) ){
  
  # ind efficacy
  p_e_i_diff <- rowDiff(get(paste0("samp_bb_e_",s_id)))
  
  # ind safety
  p_s_i_diff <- rowDiff(get(paste0("samp_bb_s_",s_id)))
  
  # joint mod
  samp_bb_jnt <- get(paste0("samp_bb_jnt_",s_id))
  
  # joint mod efficacy
  p_e_j_diff <-rowDiff(samp_bb_jnt[,c("p_e[1]","p_e[2]")])
  
  # joint mod safety
  p_s_j_diff <- rowDiff(samp_bb_jnt[,c("p_s[1]","p_s[2]")])
  
  #! order matters here because of labels
  diffs <- c("p_e_i_diff", "p_e_j_diff", "p_s_i_diff", "p_s_j_diff")
  diffs_ci <- lapply(diffs, function(x) quantile(get(x), qtile))
  diffs_ci_vec <- unlist(diffs_ci) 
  names(diffs_ci_vec) <- paste0(sort(rep(diffs,length(qtile))),"_q",qtile)
  
  return( c(s_id=s_id, diffs_ci_vec) )
}

# function to reshape combined output of ci_bb_sim
reshaper1 <- function(ds){
mlt0 <- reshape2::melt(ds, id.vars="s_id")
vars0 <- reshape2::colsplit(mlt0$variable,"_",names=c("p","outcome","mod","diff","qtile"))
mlt1 <- cbind(mlt0,vars0)
return(reshape2::dcast(mlt1, s_id+outcome+mod~qtile))
}


# function to get posterior probabilities of marginal events
pp_bb_sim <- function(s_id, grd = seq(0,1,by=0.01) ){
  
  # ind efficacy
  p_e_i_diff <- rowDiff(get(paste0("samp_bb_e_",s_id)))
  pp_e_i_gt <- sapply(grd,function(g) mean(p_e_i_diff > g)) # posterior prob efficacy diff greater than g
  
  # ind safety
  p_s_i_diff <- rowDiff(get(paste0("samp_bb_s_",s_id)))
  pp_s_i_gt <- sapply(grd,function(g) mean(p_s_i_diff > g))
  
  # joint mod
  samp_bb_jnt <- get(paste0("samp_bb_jnt_",s_id))
  
  # joint mod efficacy
  p_e_j_diff <-rowDiff(samp_bb_jnt[,c("p_e[1]","p_e[2]")])
  pp_e_j_gt <- sapply(grd,function(g) mean(p_e_j_diff > g))
  
  # joint mod safety
  p_s_j_diff <- rowDiff(samp_bb_jnt[,c("p_s[1]","p_s[2]")])
  pp_s_j_gt <- sapply(grd,function(g) mean(p_s_j_diff > g))

  #! order matters here because of labels
  postps <- c("pp_e_i_gt", "pp_e_j_gt", "pp_s_i_gt", "pp_s_j_gt")
  postps_ls <- lapply(postps, function(x) get(x))
  postps_vec <- unlist(postps_ls) 
  
  names(postps_vec) <- paste0(sort(rep(postps,length(grd))),"_g_",grd)
  
  c(s_id=s_id,postps_vec)
}

# function to reshape combined output of pp_bb_sim
reshaper2 <- function(ds){
mlt0 <- reshape2::melt(ds, id.vars="s_id", value.name="pp_gt_g")
vars0 <- reshape2::colsplit(mlt0$variable,"_",names=c("pp","outcome","mod","foo","bar","g"))
mlt1 <- cbind(mlt0,vars0)
return(mlt1[,c("s_id","outcome","mod","g","pp_gt_g")])
}



# function to get posterior probabilities of joint outcome
pp_jnt_bb_sim <- function(s_id, grd = seq(0,1,by=0.01) ){
  
  # ind efficacy
  p_e_i_diff <- rowDiff(get(paste0("samp_bb_e_",s_id)))
  
  # ind safety
  p_s_i_diff <- rowDiff(get(paste0("samp_bb_s_",s_id)))
  
  pp_jnt_i_gt <- sapply(grd,function(g) mean(p_e_i_diff > g & p_s_i_diff > g)) 
  
  # joint mod
  samp_bb_jnt <- get(paste0("samp_bb_jnt_",s_id))
  
  # joint mod efficacy
  p_e_j_diff <-rowDiff(samp_bb_jnt[,c("p_e[1]","p_e[2]")])
  
  # joint mod safety
  p_s_j_diff <- rowDiff(samp_bb_jnt[,c("p_s[1]","p_s[2]")])
  
  pp_jnt_j_gt <- sapply(grd,function(g) mean(p_s_j_diff > g & p_e_j_diff > g))

  
  #! order matters here because of labels
  postps <- c("pp_jnt_i_gt", "pp_jnt_j_gt")
  postps_ls <- lapply(postps , function(x) get(x))
  postps_vec <- unlist(postps_ls) 
  
  names(postps_vec) <- paste0(sort(rep(postps, length(grd))),"_g_",grd)
  
  c(s_id=s_id,postps_vec)
}

# function to reshape combined output of pp_jnt_bb_sim
reshaper3 <- function(ds){
mlt0 <- reshape2::melt(ds, id.vars="s_id", value.name="pp_gt_g")
vars0 <- reshape2::colsplit(mlt0$variable,"_",names=c("pp","jnt","mod","foo","bar","g"))
mlt1 <- cbind(mlt0,vars0)
return(mlt1[,c("s_id","jnt","mod","g","pp_gt_g")])
}


# function to get correlation posterior eff risk diff and safety risk diff
corr_bb_sim <- function(s_id){
  
  # ind efficacy
  p_e_i_diff <- rowDiff(get(paste0("samp_bb_e_",s_id)))
  
  # ind safety
  p_s_i_diff <- rowDiff(get(paste0("samp_bb_s_",s_id)))
  
  corr_i <- cor(p_e_i_diff, p_s_i_diff) 
  
  # joint mod
  samp_bb_jnt <- get(paste0("samp_bb_jnt_",s_id))
  
  # joint mod efficacy
  p_e_j_diff <-rowDiff(samp_bb_jnt[,c("p_e[1]","p_e[2]")])
  
  # joint mod safety
  p_s_j_diff <- rowDiff(samp_bb_jnt[,c("p_s[1]","p_s[2]")])
  
  corr_j <- cor(p_e_j_diff, p_s_j_diff) 
  
  c(s_id=s_id,corr_ind = corr_i,corr_jnt=corr_j)
}


# function to reshape combined output of pp_jnt_bb_sim
reshaper4 <- function(ds){
mlt0 <- reshape2::melt(ds, id.vars="s_id", value.name="corr")
vars0 <- reshape2::colsplit(mlt0$variable,"_",names=c("corr","mod"))
mlt1 <- cbind(mlt0,vars0)
return(mlt1[,c("s_id","mod","corr")])
}

# function to batch read-in and processing of data
# batches are deleted after loading
doBatch <- function(ntotal, batchsize){
  batches<-ntotal/batchsize
  
  # get batch sequence
  ss0<-c(seq(1,ntotal,by=batchsize),ntotal-1)
  ss2 <- ss0[-1]+1
  ss1 <- rev(rev(ss0)[-1])
  fn<-c(ss1[1], rev(rev(ss2)[-1]))
  sn<-c(ss1[-1], rev(ss2)[1])
  
  for (i in 1:batches){
    
    for (s_id in fn[i]:sn[i]){
      load(file.path(wdir,"sims","bb",paste0("bb_sim_",s_id,".RData")), .GlobalEnv)
    }
    
    # ci_bb_sim for batch - credible intervals
    assign(paste0("bb_ci_dat0_",i),
             data.frame(do.call(rbind, lapply(fn[i]:sn[i], ci_bb_sim))),
             envir=.GlobalEnv)
    
    # pp_bb_sim for batch - posterior prob gt grid
    assign(paste0("bb_pp_dat0_",i),
             data.frame(do.call(rbind, lapply(fn[i]:sn[i], pp_bb_sim))),
             envir=.GlobalEnv)
    
    # pp_jnt_bb_sim for batch - posterior prob gt grid
    assign(paste0("bb_pp_jnt_dat0_",i),
             data.frame(do.call(rbind, lapply(fn[i]:sn[i], pp_jnt_bb_sim))),
             envir=.GlobalEnv)
    
    assign(paste0("corr_dat0_",i),
           data.frame(do.call(rbind, lapply(fn[i]:sn[i], corr_bb_sim))),
           envir=.GlobalEnv)
    
    
    # remove loaded datasets, samples, summaries and sim params for batch
    rm(list=ls(envir=.GlobalEnv)[grep("dat_bb_|samp_bb_|summ_bb_|sim_params_",
                                      ls(envir=.GlobalEnv))], envir=.GlobalEnv)
  }
  
}

make function looking at correlation between p_e[2]- p_e[1] and p_s[2] - p_s[1]

# old code to make pp_jnt_bb_sim
for (i in 1:5){
  load(file.path(wdir,"sims","bb",paste0("bb_sim_",i,".RData")))
}

# function to get correlation posterior eff risk diff and safety risk diff
corr_bb_sim <- function(s_id){
  
  # ind efficacy
  p_e_i_diff <- rowDiff(get(paste0("samp_bb_e_",s_id)))
  
  # ind safety
  p_s_i_diff <- rowDiff(get(paste0("samp_bb_s_",s_id)))
  
  corr_i <- cor(p_e_i_diff, p_s_i_diff) 
  
  # joint mod
  samp_bb_jnt <- get(paste0("samp_bb_jnt_",s_id))
  
  # joint mod efficacy
  p_e_j_diff <-rowDiff(samp_bb_jnt[,c("p_e[1]","p_e[2]")])
  
  # joint mod safety
  p_s_j_diff <- rowDiff(samp_bb_jnt[,c("p_s[1]","p_s[2]")])
  
  corr_j <- cor(p_e_j_diff, p_s_j_diff) 
  
  c(s_id=s_id,corr_ind = corr_i,corr_jnt=corr_j)
}

corr_bb_sim(1)
corr_bb_sim(2)
corr_bb_sim(3)
corr_bb_sim(4)


corr_dat0<- data.frame(do.call(rbind, lapply(1:5, corr_bb_sim)))

ds<-corr_dat0


# function to reshape combined output of pp_jnt_bb_sim
reshaper4 <- function(ds){
mlt0 <- reshape2::melt(ds, id.vars="s_id", value.name="corr")
vars0 <- reshape2::colsplit(mlt0$variable,"_",names=c("corr","mod"))
mlt1 <- cbind(mlt0,vars0)
return(mlt1[,c("s_id","mod","corr")])
}

merge(bb_sim_array_ord_run, corr_dat0, by.x="sim_id",by.y="s_id")
# old code to make pp_jnt_bb_sim
for (i in 1:5){
  load(file.path(wdir,"sims","bb",paste0("bb_sim_",i,".RData")))
}

# function to get posterior probabilities of joint outcome
pp_jnt_bb_sim <- function(s_id, grd = seq(0,1,by=0.01) ){
  
  # ind efficacy
  p_e_i_diff <- rowDiff(get(paste0("samp_bb_e_",s_id)))
  
  # ind safety
  p_s_i_diff <- rowDiff(get(paste0("samp_bb_s_",s_id)))
  
  pp_jnt_i_gt <- sapply(grd,function(g) mean(p_e_i_diff > g & p_s_i_diff > g)) 
  
  # joint mod
  samp_bb_jnt <- get(paste0("samp_bb_jnt_",s_id))
  
  # joint mod efficacy
  p_e_j_diff <-rowDiff(samp_bb_jnt[,c("p_e[1]","p_e[2]")])
  
  # joint mod safety
  p_s_j_diff <- rowDiff(samp_bb_jnt[,c("p_s[1]","p_s[2]")])
  
  pp_jnt_j_gt <- sapply(grd,function(g) mean(p_s_j_diff > g & p_e_j_diff > g))

  
  #! order matters here because of labels
  postps <- c("pp_jnt_i_gt", "pp_jnt_j_gt")
  postps_ls <- lapply(postps , function(x) get(x))
  postps_vec <- unlist(postps_ls) 
  
  names(postps_vec) <- paste0(sort(rep(postps, length(grd))),"_g_",grd)
  
  c(s_id=s_id,postps_vec)
}

ds<-pp_jnt_bb_sim(1)

bb_pp_jnt_dat0<- data.frame(do.call(rbind, lapply(1:5, pp_jnt_bb_sim)))

ds <- bb_pp_jnt_dat0
mlt0 <- reshape2::melt(ds, id.vars="s_id", value.name="pp_gt_g")
vars0 <- reshape2::colsplit(mlt0$variable,"_",names=c("pp","jnt","mod","foo","bar","g"))
mlt1 <- cbind(mlt0,vars0)

# function to reshape combined output of pp_bb_sim
reshaper3 <- function(ds){
mlt0 <- reshape2::melt(ds, id.vars="s_id", value.name="pp_gt_g")
vars0 <- reshape2::colsplit(mlt0$variable,"_",names=c("pp","jnt","mod","foo","bar","g"))
mlt1 <- cbind(mlt0,vars0)
return(mlt1[,c("s_id","jnt","mod","g","pp_gt_g")])
}

reshaper3(ds)
# read in array with sim params
bb_sim_array <- readRDS(file.path(wdir,"bb_sim_array.rds"))
bb_sim_array_ord<-bb_sim_array[order(bb_sim_array$n,bb_sim_array$rho_2,bb_sim_array$p_e2,bb_sim_array$p_s2,bb_sim_array$rep_id),]

nreps<-100
nscn<-nrow(bb_sim_array)/nreps
sc_id<-sort(rep(1:nscn,nreps))
bb_sim_array_ord$scn_id <- sc_id

# keep only ones where sim has been run
nsims<-2000
bb_sim_array_ord_run <- subset(bb_sim_array_ord, sim_id<=nsims)
rm(bb_sim_array, bb_sim_array_ord, nreps, nscn, sc_id)

# placebo group
p_e1 <- 0.2 # prob effective
p_s1 <- 0.1 # prob AE
rho_1 <- 0.1 # tetrachoric correlation, can estimate with polycor::polychor
# load and process sims in batches
bsize<-200

#doBatch(100, 20)

doBatch(nsims, bsize)

bb_ci_vec<-ls()[grep("bb_ci_dat0_",ls())]
bb_ci_lst<-sapply(bb_ci_vec, get,simplify = FALSE)
bb_ci_dat0<-do.call(rbind,bb_ci_lst)
rm(list=bb_ci_vec)
rm(bb_ci_lst,bb_ci_vec)

bb_ci_dat1 <- reshaper1(bb_ci_dat0)
bb_ci_dat<-merge(bb_sim_array_ord_run, bb_ci_dat1, by.x="sim_id",by.y="s_id")
rm(bb_ci_dat0,bb_ci_dat1)


bb_pp_vec<-ls()[grep("bb_pp_dat0_",ls())]
bb_pp_lst<-sapply(bb_pp_vec, get,simplify = FALSE)
bb_pp_dat0<-do.call(rbind,bb_pp_lst)
rm(list=bb_pp_vec)
rm(bb_pp_lst,bb_pp_vec)

bb_pp_dat1 <- reshaper2(bb_pp_dat0)
bb_pp_dat<-merge(bb_sim_array_ord_run, bb_pp_dat1, by.x="sim_id",by.y="s_id")
rm(bb_pp_dat0,bb_pp_dat1)


bb_pp_jnt_vec<-ls()[grep("bb_pp_jnt_dat0_",ls())]
bb_pp_jnt_lst<-sapply(bb_pp_jnt_vec, get,simplify = FALSE)
bb_pp_jnt_dat0<-do.call(rbind,bb_pp_jnt_lst)
rm(list=bb_pp_jnt_vec)
rm(bb_pp_jnt_lst,bb_pp_jnt_vec)

bb_pp_jnt_dat1 <- reshaper3(bb_pp_jnt_dat0)
bb_pp_jnt_dat<-merge(bb_sim_array_ord_run, bb_pp_jnt_dat1, by.x="sim_id",by.y="s_id")
rm(bb_pp_jnt_dat0,bb_pp_jnt_dat1)

corr_vec<-ls()[grep("corr_dat0_",ls())]
corr_lst<-sapply(corr_vec, get,simplify = FALSE)
corr_dat0<-do.call(rbind,corr_lst)
rm(list=corr_vec)
rm(corr_lst,corr_vec)

corr_dat1 <- reshaper4(corr_dat0)
corr_dat<-merge(bb_sim_array_ord_run, corr_dat1, by.x="sim_id",by.y="s_id")
rm(corr_dat0,corr_dat1)

plots of credible intervals

library(ggstance)
## 
## Attaching package: 'ggstance'
## The following objects are masked from 'package:ggplot2':
## 
##     geom_errorbarh, GeomErrorbarh
if (0){
#too much!
ggplot(data=bb_ci_dat, aes(x=q0.5, y=outcome, group=mod)) + geom_point(aes(shape=mod), position=position_dodgev(height=0.2)) + facet_grid(n~rho_2+p_e2+p_s2)

ggplot(data=bb_ci_dat,aes(x=q0.5, xmin=q0.05, xmax=q0.95, y=outcome, color=mod)) + 
  geom_point(aes(shape=mod), position=position_dodgev(height=0.4), alpha=0.4) +
  geom_linerangeh(position=position_dodgev(height=0.4), alpha=0.4) + facet_grid(n~rho_2+p_e2+p_s2)

# better
ggplot(data=subset(bb_ci_dat,rho_2==0.35), aes(x=q0.5, y=outcome, group=mod)) + geom_point(aes(shape=mod), position=position_dodgev(height=0.2)) + facet_grid(n~p_e2+p_s2)

ggplot(data=subset(bb_ci_dat,rho_2==0.35),aes(x=q0.5, xmin=q0.05, xmax=q0.95, y=outcome, color=mod)) + geom_point(aes(shape=mod), position=position_dodgev(height=0.4), alpha=0.4) +
  geom_linerangeh(position=position_dodgev(height=0.4), alpha=0.4) + facet_grid(n~p_e2+p_s2)
}

if (1){
ggplot(data=subset(bb_ci_dat,rho_2==0.35 & n==50),aes(x=q0.5, xmin=q0.05, xmax=q0.95, y=outcome, color=mod)) +   geom_point(aes(shape=mod), position=position_dodgev(height=0.4), alpha=0.1) +
  geom_linerangeh(position=position_dodgev(height=0.4), alpha=0.25) + 
    facet_grid(p_e2~p_s2,labeller = "label_both")

ggplot(data=subset(bb_ci_dat,rho_2==0.35 & n==100),aes(x=q0.5, xmin=q0.05, xmax=q0.95, y=outcome, color=mod)) +   geom_point(aes(shape=mod), position=position_dodgev(height=0.4), alpha=0.1) +
  geom_linerangeh(position=position_dodgev(height=0.4), alpha=0.25) + facet_grid(p_e2~p_s2,labeller = "label_both")

ggplot(data=subset(bb_ci_dat,rho_2==0.35 & n==200),aes(x=q0.5, xmin=q0.05, xmax=q0.95, y=outcome, color=mod)) +   geom_point(aes(shape=mod), position=position_dodgev(height=0.4), alpha=0.1) +
  geom_linerangeh(position=position_dodgev(height=0.4), alpha=0.25) + facet_grid(p_e2~p_s2,labeller = "label_both")

# with 400 really tight

}

if (1){
  ggplot(data=subset(bb_ci_dat,rho_2==0.1 & n==100),aes(x=q0.5, xmin=q0.05, xmax=q0.95, y=outcome, color=mod)) +   geom_point(aes(shape=mod), position=position_dodgev(height=0.4), alpha=0.1) +
  geom_linerangeh(position=position_dodgev(height=0.4), alpha=0.25) + facet_grid(p_e2~p_s2,labeller = "label_both")
  
  ggplot(data=subset(bb_ci_dat,rho_2==0.35 & n==100),aes(x=q0.5, xmin=q0.05, xmax=q0.95, y=outcome, color=mod)) +   geom_point(aes(shape=mod), position=position_dodgev(height=0.4), alpha=0.1) +
  geom_linerangeh(position=position_dodgev(height=0.4), alpha=0.25) + facet_grid(p_e2~p_s2,labeller = "label_both")

    ggplot(data=subset(bb_ci_dat,rho_2==0.6 & n==100),aes(x=q0.5, xmin=q0.05, xmax=q0.95, y=outcome, color=mod)) +   geom_point(aes(shape=mod), position=position_dodgev(height=0.4), alpha=0.1) +  geom_linerangeh(position=position_dodgev(height=0.4), alpha=0.25) + facet_grid(p_e2~p_s2,labeller = "label_both")
  }

# test with single scenerio
subset(bb_ci_dat,scn_id==72)

ggplot(data=subset(bb_ci_dat,scn_id==72), aes(x=q0.5, y=outcome, group=mod)) + geom_point(aes(shape=mod), position=position_dodgev(height=0.2))

ggplot(data=subset(bb_ci_dat,scn_id==72),aes(x=q0.5, xmin=q0.05, xmax=q0.95, y=outcome, color=mod)) + 
  geom_point(aes(shape=mod), position=position_dodgev(height=0.4), alpha=0.4) +
  geom_linerangeh(position=position_dodgev(height=0.4), alpha=0.4) 

#+    geom_point(data=data.frame(x=c(p_e1, p_e2, p_s1, p_s2),y=1:4), aes(x=x,y=y),
#             color="blue", size=2, inherit.aes = FALSE)

plots of ind posterior prob greater than g

bb_pp_dat_srt<-bb_pp_dat %>% arrange(rep_id,outcome,mod,g)
ss<-subset(bb_pp_dat,scn_id==72)
table(ss$mod,ss$outcome,ss$rep_id)
## , ,  = 5
## 
##    
##       e   s
##   i 101 101
##   j 101 101
## 
## , ,  = 6
## 
##    
##       e   s
##   i 101 101
##   j 101 101
## 
## , ,  = 16
## 
##    
##       e   s
##   i 101 101
##   j 101 101
## 
## , ,  = 18
## 
##    
##       e   s
##   i 101 101
##   j 101 101
## 
## , ,  = 27
## 
##    
##       e   s
##   i 101 101
##   j 101 101
## 
## , ,  = 32
## 
##    
##       e   s
##   i 101 101
##   j 101 101
## 
## , ,  = 39
## 
##    
##       e   s
##   i 101 101
##   j 101 101
## 
## , ,  = 52
## 
##    
##       e   s
##   i 101 101
##   j 101 101
## 
## , ,  = 55
## 
##    
##       e   s
##   i 101 101
##   j 101 101
## 
## , ,  = 63
## 
##    
##       e   s
##   i 101 101
##   j 101 101
## 
## , ,  = 72
## 
##    
##       e   s
##   i 101 101
##   j 101 101
## 
## , ,  = 73
## 
##    
##       e   s
##   i 101 101
##   j 101 101
## 
## , ,  = 74
## 
##    
##       e   s
##   i 101 101
##   j 101 101
## 
## , ,  = 75
## 
##    
##       e   s
##   i 101 101
##   j 101 101
## 
## , ,  = 77
## 
##    
##       e   s
##   i 101 101
##   j 101 101
## 
## , ,  = 87
## 
##    
##       e   s
##   i 101 101
##   j 101 101
## 
## , ,  = 88
## 
##    
##       e   s
##   i 101 101
##   j 101 101
## 
## , ,  = 91
## 
##    
##       e   s
##   i 101 101
##   j 101 101
## 
## , ,  = 94
## 
##    
##       e   s
##   i 101 101
##   j 101 101
## 
## , ,  = 97
## 
##    
##       e   s
##   i 101 101
##   j 101 101
## 
## , ,  = 99
## 
##    
##       e   s
##   i 101 101
##   j 101 101
ggplot(data=subset(bb_pp_dat_srt,scn_id==72 & rep_id==99), aes(x=g, y=pp_gt_g, linetype=outcome)) + geom_line(aes(col=mod))

ggplot(data=subset(bb_pp_dat_srt,scn_id==72 & rep_id==97), aes(x=g, y=pp_gt_g, linetype=outcome)) + geom_line(aes(col=mod))

ggplot(data=subset(bb_pp_dat_srt,scn_id==72), aes(x=g, y=pp_gt_g, group=rep_id)) + geom_line(alpha=0.5) + facet_grid(mod~outcome)

# 
ggplot(data=subset(bb_pp_dat_srt,scn_id==72), 
       aes(x=g, y=pp_gt_g, color=mod, group=interaction(mod,rep_id))) + 
  geom_line(alpha=0.5) + 
  facet_grid(outcome~.)

ggplot(data=subset(bb_pp_dat,rho_2==0.35 & n==50 & outcome=="e"), 
       aes(x=g, y=pp_gt_g, color=mod, group=interaction(mod,rep_id))) + 
  geom_line(alpha=0.5) + facet_grid(p_e2~p_s2,labeller = "label_both")

ggplot(data=subset(bb_pp_dat,rho_2==0.35 & n==50), aes(x=g, y=pp_gt_g, group=rep_id)) + geom_line(aes(col=mod)) + facet_grid(outcome+p_e2~p_s2,labeller = "label_both")

plots of joint posterior prob greater than g

# 
ggplot(data=subset(bb_pp_jnt_dat,scn_id==72), 
       aes(x=g, y=pp_gt_g, color=mod, group=interaction(mod,rep_id))) + 
  geom_line(alpha=0.5) 

ggplot(data=subset(bb_pp_jnt_dat,scn_id==1), 
       aes(x=g, y=pp_gt_g, color=mod, group=interaction(mod,rep_id))) + 
  geom_line(alpha=0.5) 

ggplot(data=subset(bb_pp_jnt_dat, rho_2==0.1 & n==100), 
       aes(x=g, y=pp_gt_g, color=mod, group=interaction(mod,rep_id))) + 
  geom_line(alpha=0.5) + facet_grid(p_e2~p_s2,labeller = "label_both")

ggplot(data=subset(bb_pp_jnt_dat, rho_2==0.35 & n==100), 
       aes(x=g, y=pp_gt_g, color=mod, group=interaction(mod,rep_id))) + 
  geom_line(alpha=0.5) + facet_grid(p_e2~p_s2,labeller = "label_both")

ggplot(data=subset(bb_pp_jnt_dat, rho_2==0.6 & n==100), 
       aes(x=g, y=pp_gt_g, color=mod, group=interaction(mod,rep_id))) + 
  geom_line(alpha=0.5) + facet_grid(p_e2~p_s2,labeller = "label_both")

plots of corr

ggplot(data=subset(corr_dat,scn_id==72), aes(x=corr, fill=mod)) + geom_density()

ggplot(data=subset(corr_dat,rho_2==0.35 & n==100),aes(x=corr, fill=mod)) +  geom_density( alpha=0.1) + facet_grid(p_e2~p_s2,labeller = "label_both")

ggplot(data=subset(corr_dat,p_e2==0.2 & p_s2==0.1),aes(x=corr, y=mod)) +  geom_violinh(alpha=0.1) + facet_grid(n~rho_2,labeller = "label_both")

ggplot(data=subset(corr_dat,p_e2==0.2 & p_s2==0.1),aes(x=corr, fill=mod)) +  geom_density( alpha=0.1) + facet_grid(n~rho_2,labeller = "label_both")

ggplot(data=subset(corr_dat,p_e2==0.2 & p_s2==0.4),aes(x=corr, fill=mod)) +  geom_density( alpha=0.1) + facet_grid(n~rho_2,labeller = "label_both")

ggplot(data=subset(corr_dat,p_e2==0.2 & p_s2==0.7),aes(x=corr, fill=mod)) +  geom_density( alpha=0.1) + facet_grid(n~rho_2,labeller = "label_both")

ggplot(data=subset(corr_dat,p_e2==0.5 & p_s2==0.1),aes(x=corr, fill=mod)) +  geom_density( alpha=0.1) + facet_grid(n~rho_2,labeller = "label_both")

ggplot(data=subset(corr_dat,p_e2==0.5 & p_s2==0.4),aes(x=corr, fill=mod)) +  geom_density( alpha=0.1) + facet_grid(n~rho_2,labeller = "label_both")

ggplot(data=subset(corr_dat,p_e2==0.5 & p_s2==0.7),aes(x=corr, fill=mod)) +  geom_density( alpha=0.1) + facet_grid(n~rho_2,labeller = "label_both")

ggplot(data=subset(corr_dat,p_e2==0.8 & p_s2==0.1),aes(x=corr, fill=mod)) +  geom_density( alpha=0.1) + facet_grid(n~rho_2,labeller = "label_both")

ggplot(data=subset(corr_dat,p_e2==0.8 & p_s2==0.4),aes(x=corr, fill=mod)) +  geom_density( alpha=0.1) + facet_grid(n~rho_2,labeller = "label_both")

ggplot(data=subset(corr_dat,p_e2==0.8 & p_s2==0.7),aes(x=corr, fill=mod)) +  geom_density( alpha=0.1) + facet_grid(n~rho_2,labeller = "label_both")